%% 初始化设置
% 重置路径避免冲突
restoredefaultpath;

% 添加工具箱路径（使用安全方式）
ftPath = 'C:\Users\custs\workspace\matlab\congcong\lib\fieldtrip-master';
addpath(ftPath);
ft_defaults;  % FieldTrip初始化

% 设置数据路径
dataDir = 'C:\Users\custs\workspace\matlab\congcong\data\aMCI\clean';
outputDir = 'C:\Users\custs\workspace\matlab\congcong\data\mat';
layoutFile = 'C:\Users\custs\workspace\matlab\congcong\layout\biosemi64.lay';

% 确保输出目录存在
if ~exist(outputDir, 'dir')
    mkdir(outputDir);
end

% 被试列表
subjectIDs = {'103', '106', '204'};
Ns = length(subjectIDs);

%% 主分析循环
for subIdx = 1:Ns
    fprintf('\n===== 处理被试 %d/%d: %s =====\n', subIdx, Ns, subjectIDs{subIdx});
    
    try
        %% 加载数据
        load(fullfile(dataDir, [subjectIDs{subIdx} '_amciCleanData.mat']));
        
        %% 时频分析配置
        cfg = struct();
        cfg.output = 'fourier';
        cfg.pad = 'nextpow2';
        cfg.channel = 'all';
        cfg.method = 'mtmconvol';
        cfg.taper = 'hanning';
        cfg.foi = 1:0.05:30;  % 频率范围:1-30Hz,步长0.05Hz
        cfg.t_ftimwin = 0.5 * ones(size(cfg.foi)); % 时间窗固定为0.5秒
        cfg.toi = -0.7:0.001:1.5; % 时间点:-0.7到1.5秒,步长1ms
        
        %% 计算时频表征
        fprintf('计算条件A的时频表征...\n');
        fourier_A = ft_freqanalysis(cfg, condA);
        
        fprintf('计算条件B的时频表征...\n');
        fourier_B = ft_freqanalysis(cfg, condB);
        
        %% 计算ITC(跨试验相位一致性)
        fprintf('计算ITC...\n');
        
        % 条件A
        itc_A = struct();
        itc_A.label = fourier_A.label;
        itc_A.freq = fourier_A.freq;
        itc_A.time = fourier_A.time;
        itc_A.dimord = 'chan_freq_time';
        F_A = fourier_A.fourierspctrm;
        itc_A.itc = squeeze(abs(mean(F_A./abs(F_A), 1)));
        
        % 条件B
        itc_B = itc_A;
        F_B = fourier_B.fourierspctrm;
        itc_B.itc = squeeze(abs(mean(F_B./abs(F_B), 1)));
        
        %% 保存结果
        save(fullfile(outputDir, [subjectIDs{subIdx} '.mat']), 'itc_A', 'itc_B');
        fprintf('结果已保存\n');
        
    catch ME
        fprintf(2, '处理被试 %s 时出错: %s\n', subjectIDs{subIdx}, ME.message);
        fprintf(2, '发生在: %s (行 %d)\n', ME.stack(1).name, ME.stack(1).line);
        continue;
    end
end

%% 组水平分析
fprintf('\n===== 开始组水平分析 =====\n');

% 准备数据结构
allsub_A = cell(1, Ns);
allsub_B = cell(1, Ns);

for subIdx = 1:Ns
    load(fullfile(outputDir, [subjectIDs{subIdx} '.mat']));
    allsub_A{subIdx} = itc_A;
    allsub_B{subIdx} = itc_B;
end

% 计算组平均
cfg = struct();
cfg.parameter = 'itc';
cfg.channel = 'all';
cfg.toilim = 'all';
cfg.foilim = 'all';

fprintf('计算组平均ITC...\n');
itc_A_grand = ft_freqgrandaverage(cfg, allsub_A{:});
itc_B_grand = ft_freqgrandaverage(cfg, allsub_B{:});

% 条件间差异
A_vs_B = itc_A_grand;
A_vs_B.itc = itc_A_grand.itc - itc_B_grand.itc;

%% 可视化结果
% 配置通用绘图参数
plot_cfg = struct();
plot_cfg.parameter = 'itc';
plot_cfg.showlabels = 'yes';
plot_cfg.layout = layoutFile;
plot_cfg.colorbar = 'yes';

% 1. 组平均时频图
figure('Name', '组平均ITC - 条件A', 'Position', [100 100 800 600]);
ft_multiplotTFR(plot_cfg, itc_A_grand);
title('组平均ITC - 条件A');

figure('Name', '组平均ITC - 条件B', 'Position', [100 100 800 600]);
ft_multiplotTFR(plot_cfg, itc_B_grand);
title('组平均ITC - 条件B');

% 2. 差异地形图
topo_cfg = plot_cfg;
topo_cfg.comment = 'no';
topo_cfg.xlim = [0 1.5]; % 时间范围
topo_cfg.ylim = [8 13];  % Alpha频段
topo_cfg.zlim = 'maxabs'; % 自动设置对称的颜色范围

figure('Name', 'ITC差异地形图(Alpha波段)', 'Position', [100 100 800 600]);
ft_topoplotTFR(topo_cfg, A_vs_B);
title('条件A - 条件B (Alpha波段: 8-13Hz)');

%% ROI分析
roi_channels = {'P3', 'P4', 'P5', 'P6', 'PO3', 'PO4', 'PO7', 'PO8'};
timewin = [0.15 1.5]; % 分析的时间窗口
freqwin = [8 13];    % Alpha频段

% 预分配结果矩阵
power_results = zeros(Ns, 2); % 每个被试两个条件的功率

for ch = 1:length(roi_channels)
    chan = roi_channels{ch};
    
    % 检查电极是否存在
    if ~any(strcmpi(chan, itc_A_grand.label))
        fprintf('警告: 电极 %s 不存在于数据中\n', chan);
        continue;
    end
    
    % 创建单电极图
    fig = figure('Name', ['电极 ' chan ' 分析'], 'Position', [100 100 1000 800]);
    
    % 1. 条件A
    cfg = plot_cfg;
    cfg.figure = subplot(2,2,1);
    cfg.channel = chan;
    ft_singleplotTFR(cfg, itc_A_grand);
    title(['条件A - ' chan]);
    
    % 2. 条件B
    cfg.figure = subplot(2,2,2);
    ft_singleplotTFR(cfg, itc_B_grand);
    title(['条件B - ' chan]);
    
    % 3. 差异
    cfg.figure = subplot(2,2,3);
    ft_singleplotTFR(cfg, A_vs_B);
    title(['差异(A-B) - ' chan]);
    
    % 4. 统计分析
    cfg.figure = subplot(2,2,4);
    
    % 提取ROI数据
    time_idx = dsearchn(itc_A_grand.time', timewin');
    freq_idx = dsearchn(itc_A_grand.freq', freqwin');
    chan_idx = find(strcmpi(chan, itc_A_grand.label));
    
    % 计算每个被试在ROI内的平均ITC
    for subIdx = 1:Ns
        powA = allsub_A{subIdx}.itc(chan_idx, freq_idx(1):freq_idx(2), time_idx(1):time_idx(2));
        powB = allsub_B{subIdx}.itc(chan_idx, freq_idx(1):freq_idx(2), time_idx(1):time_idx(2));
        
        power_results(subIdx, 1) = mean(powA(:));
        power_results(subIdx, 2) = mean(powB(:));
    end
    
    % 配对t检验
    [h, p, ci, stats] = ttest(power_results(:,1), power_results(:,2));
    
    % 绘制柱状图
    bar(mean(power_results));
    hold on;
    errorbar(1:2, mean(power_results), std(power_results)/sqrt(Ns), 'k.', 'LineWidth', 1.5);
    title(sprintf('%s: p = %.4f (t = %.2f)', chan, p, stats.tstat));
    xticks(1:2);
    xticklabels({'条件A', '条件B'});
    ylabel('平均ITC');
    
    % 保存图形
    saveas(fig, fullfile(outputDir, [chan '_analysis.png']));
end

%% 导出数据用于进一步分析
% 选择中央区域电极
central_channels = {'FZ', 'FCZ', 'CZ'};
chan_idx = find(ismember(itc_A_grand.label, central_channels));

% 定义分析窗口
timewin = [-0.4 1.2];
time_idx = dsearchn(itc_A_grand.time', timewin');
freqwin = [1 30];  % 全频段
freq_idx = dsearchn(itc_A_grand.freq', freqwin');

% 提取数据
results = zeros(Ns, 2); % 每个被试两个条件

for subIdx = 1:Ns
    powA = allsub_A{subIdx}.itc(chan_idx, freq_idx(1):freq_idx(2), time_idx(1):time_idx(2));
    powB = allsub_B{subIdx}.itc(chan_idx, freq_idx(1):freq_idx(2), time_idx(1):time_idx(2));
    
    results(subIdx, 1) = mean(powA(:));
    results(subIdx, 2) = mean(powB(:));
end

% 保存为CSV文件
csvwrite(fullfile(outputDir, 'itc_results.csv'), results);
fprintf('\n===== 分析完成 =====\n');